package ssau.knyazev.fft.modules.real;

import ssau.knyazev.modules.Generator;

public class RealCFFTModule {

	protected RealCFFTModule(){ }
	

	public static float[][] _fft(float[][] src){
		int N = src.length;
		int size = 0;
		for (int x = src.length; x > 1; x = x >> 1){
			size++;
		}
		inverse(src);
		src = transformToGyperComplexMatrix(src);
		float[] wi_list = new float[src.length*4];
		float[] wj_list = new float[src.length*4];
		for (int i = 0; i < src.length; i++) {
			double power = -2*Math.PI*i/N;
			wi_list[4*i] = (float)Math.cos(power);
			wi_list[4*i+1] = (float)Math.sin(power);
			wi_list[4*i+2] = 0;
			wi_list[4*i+3] = 0;
			wj_list[4*i] = (float)Math.cos(power);
			wj_list[4*i+1] = 0;
			wj_list[4*i+2] =(float)Math.sin(power);
			wj_list[4*i+3] = 0;
		}
		float[] tmp = new float[4];
		float[] tmp2 = new float[4];
		for (int q = 0; q < size; q++) {
			int step = 1 << q;
			int matrixWidth = step << 1;
			for (int down = 0; down < N - step; down++) {
				for (int right = 0; right < N - step; right++) {
					boolean flag = (down/step) % 2 + (right/step) % 2 > 0 ? false : true;
					if (flag) {
						int wiInd = (down % matrixWidth)*(N/matrixWidth)*4;
						int wjInd = (right % matrixWidth)*(N/matrixWidth)*4;
						float[] wi = {wi_list[wiInd], wi_list[wiInd+1], wi_list[wiInd+2], wi_list[wiInd+3]};
						float[] wj = {wj_list[wjInd], wj_list[wjInd+1], wj_list[wjInd+2], wj_list[wjInd+3]};
						
						int raInd = 4*right;
						int rbInd = raInd + 1;
						int rcInd = raInd + 2;
						int rdInd = raInd + 3;
						int rsaInd = 4*(right+step);
						int rsbInd = rsaInd+1;
						int rscInd = rsaInd+2;
						int rsdInd = rsaInd+3;
						
						float a00a = src[down][raInd];
						float a00b = src[down][rbInd];
						float a00c = src[down][rcInd];
						float a00d = src[down][rdInd];
						
						float a10a = src[down][raInd];
						float a10b = src[down][rbInd];
						float a10c = src[down][rcInd];
						float a10d = src[down][rdInd];
						
						float a01a = src[down][raInd];
						float a01b = src[down][rbInd];
						float a01c = src[down][rcInd];
						float a01d = src[down][rdInd];
						
						float a11a = src[down][raInd];
						float a11b = src[down][rbInd];
						float a11c = src[down][rcInd];
						float a11d = src[down][rdInd];
						
						//XXX = XXX.add(t1.multy(a10));
						tmp[0] = wi[0]*src[down+step][raInd] - wi[1]*src[down+step][rbInd] - wi[2]*src[down+step][rcInd] + wi[3]*src[down+step][rdInd];
						a00a = a00a + tmp[0];
						a10a = a10a - tmp[0];
						a01a = a01a + tmp[0];
						a11a = a11a - tmp[0];
						
						tmp[1] = wi[0]*src[down+step][rbInd] + wi[1]*src[down+step][raInd] - wi[2]*src[down+step][rdInd] + wi[3]*src[down+step][rcInd];
						a00b = a00b + tmp[1];
						a10b = a10b - tmp[1];
						a01b = a01b + tmp[1];
						a11b = a11b - tmp[1];
						
						tmp[2] = wi[0]*src[down+step][rcInd] - wi[1]*src[down+step][rdInd] + wi[2]*src[down+step][raInd] - wi[3]*src[down+step][rbInd];
						a00c = a00c + tmp[2];
						a10c = a10c - tmp[2];
						a01c = a01c + tmp[2];
						a11c = a11c - tmp[2];
						
						tmp[3] = wi[0]*src[down+step][rdInd] + wi[1]*src[down+step][rcInd] + wi[2]*src[down+step][rbInd] + wi[3]*src[down+step][raInd];							
						a00d = a00d + tmp[3];
						a10d = a10d - tmp[3];
						a01d = a01d + tmp[3];
						a11d = a11d - tmp[3];
						
						//XXX = XXX.add(t2.multy(a01));
						tmp[0] = wj[0]*src[down][rsaInd] - wj[1]*src[down][rsbInd] - wj[2]*src[down][rscInd] + wj[3]*src[down][rsdInd];
						a00a = a00a + tmp[0];
						a10a = a10a + tmp[0];
						a01a = a01a - tmp[0];
						a11a = a11a - tmp[0];
						
						tmp[1] = wj[0]*src[down][rsbInd] + wj[1]*src[down][rsaInd] - wj[2]*src[down][rsdInd] + wj[3]*src[down][rscInd];
						a00b = a00b + tmp[1];
						a10b = a10b + tmp[1];
						a01b = a01b - tmp[1];
						a11b = a11b - tmp[1];
						
						tmp[2] = wj[0]*src[down][rscInd] - wj[1]*src[down][rsdInd] + wj[2]*src[down][rsaInd] - wj[3]*src[down][rsbInd];
						a00c = a00c + tmp[2];
						a10c = a10c + tmp[2];
						a01c = a01c - tmp[2];
						a11c = a11c - tmp[2];
						
						tmp[3] = wj[0]*src[down][rsdInd] + wj[1]*src[down][rscInd] + wj[2]*src[down][rsbInd] + wj[3]*src[down][rsaInd];	
						a00d = a00d + tmp[3];
						a10d = a10d + tmp[3];
						a01d = a01d - tmp[3];
						a11d = a11d - tmp[3];
						
						//XXX = XXX.add(t1.multy(t2.multy(a11)));
						tmp[0] = wj[0]*src[down+step][rsaInd] - wj[1]*src[down+step][rsbInd] - wj[2]*src[down+step][rscInd] + wj[3]*src[down+step][rsdInd];
						tmp[1] = wj[0]*src[down+step][rsbInd] + wj[1]*src[down+step][rsaInd] - wj[2]*src[down+step][rsdInd] + wj[3]*src[down+step][rscInd];
						tmp[2] = wj[0]*src[down+step][rscInd] - wj[1]*src[down+step][rsdInd] + wj[2]*src[down+step][rsaInd] - wj[3]*src[down+step][rsbInd];
						tmp[3] = wj[0]*src[down+step][rsdInd] + wj[1]*src[down+step][rscInd] + wj[2]*src[down+step][rsbInd] + wj[3]*src[down+step][rsaInd];	
						tmp2[0] = wi[0]*tmp[0] - wi[1]*tmp[1] - wi[2]*tmp[2] + wi[3]*tmp[3];
						tmp2[1] = wi[0]*tmp[1] + wi[1]*tmp[0] - wi[2]*tmp[3] + wi[3]*tmp[2];
						tmp2[2] = wi[0]*tmp[2] - wi[1]*tmp[3] + wi[2]*tmp[0] - wi[3]*tmp[1];
						tmp2[3] = wi[0]*tmp[3] + wi[1]*tmp[2] + wi[2]*tmp[1] + wi[3]*tmp[0];	
						
						a00a = a00a + tmp2[0];
						a00b = a00b + tmp2[1];
						a00c = a00c + tmp2[2];
						a00d = a00d + tmp2[3];
						
						a10a = a10a - tmp2[0];
						a10b = a10b - tmp2[1];
						a10c = a10c - tmp2[2];
						a10d = a10d - tmp2[3];
						
						a01a = a01a - tmp2[0];
						a01b = a01b - tmp2[1];
						a01c = a01c - tmp2[2];
						a01d = a01d - tmp2[3];
						
						a11a = a11a + tmp2[0];
						a11b = a11b + tmp2[1];
						a11c = a11c + tmp2[2];
						a11d = a11d + tmp2[3];
						
						src[down][raInd] = a00a;
						src[down][rsaInd] = a01a;
						src[down+step][raInd] = a10a;
						src[down+step][rsaInd] = a11a;
						
						src[down][rbInd] = a00b;
						src[down][rsbInd] = a01b;
						src[down+step][rbInd] = a10b;
						src[down+step][rsbInd] = a11b;
						
						src[down][rcInd] = a00c;
						src[down][rscInd] = a01c;
						src[down+step][rcInd] = a10c;
						src[down+step][rscInd] = a11c;
						
						src[down][rdInd] = a00d;
						src[down][rsdInd] = a01d;
						src[down+step][rdInd] = a10d;
						src[down+step][rsdInd] = a11d;
					
					}
				}
			}
		}
		return src;
	}

	
	public static float[][] _fftMirror(float[][] src) {
		int N = src.length;
		int size = 0;
		for (int x = src.length; x > 1; x = x >> 1){
			size++;
		}
		inverse(src);
		src = transformToGyperComplexMatrix(src);
		float[] wi_list = new float[src.length*4];
		float[] wj_list = new float[src.length*4];
		for (int i = 0; i < src.length; i++) {
			double power = -2*Math.PI*i/N;
			wi_list[4*i] = (float)Math.cos(power);
			wi_list[4*i+1] = (float)Math.sin(power);
			wi_list[4*i+2] = 0;
			wi_list[4*i+3] = 0;
			wj_list[4*i] = (float)Math.cos(power);
			wj_list[4*i+1] = 0;
			wj_list[4*i+2] =(float)Math.sin(power);
			wj_list[4*i+3] = 0;
		}
		float[] tmp = new float[4];
		float[] tmp2 = new float[4];
		for (int q = 0; q < size; q++) {
			int step = 1 << q;
			int matrixWidth = step << 1;
			int part = matrixWidth/4;
			for (int down = 0; down < N - step; down++) {
				for (int right = 0; right < N - step; right++) {
					int dmatrInd = down/matrixWidth * matrixWidth;
					int rmatrInd = right/matrixWidth * matrixWidth;
					int tdown = down - dmatrInd;
					int tright = right - rmatrInd;
					boolean flag = (tdown <= part && tright <= part);
					if (flag) {
						int wiInd = (down % matrixWidth)*(N/matrixWidth)*4;
						int wjInd = (right % matrixWidth)*(N/matrixWidth)*4;
						float[] wi = {wi_list[wiInd], wi_list[wiInd+1], wi_list[wiInd+2], wi_list[wiInd+3]};
						float[] wj = {wj_list[wjInd], wj_list[wjInd+1], wj_list[wjInd+2], wj_list[wjInd+3]};
						
						int raInd = 4*right;
						int rbInd = raInd + 1;
						int rcInd = raInd + 2;
						int rdInd = raInd + 3;
						int rsaInd = 4*(right+step);
						int rsbInd = rsaInd+1;
						int rscInd = rsaInd+2;
						int rsdInd = rsaInd+3;
						
						float a00a = src[down][raInd];
						float a00b = src[down][rbInd];
						float a00c = src[down][rcInd];
						float a00d = src[down][rdInd];
						
						float a10a = src[down][raInd];
						float a10b = src[down][rbInd];
						float a10c = src[down][rcInd];
						float a10d = src[down][rdInd];
						
						float a01a = src[down][raInd];
						float a01b = src[down][rbInd];
						float a01c = src[down][rcInd];
						float a01d = src[down][rdInd];
						
						float a11a = src[down][raInd];
						float a11b = src[down][rbInd];
						float a11c = src[down][rcInd];
						float a11d = src[down][rdInd];
						
						//XXX = XXX.add(t1.multy(a10));
						tmp[0] = wi[0]*src[down+step][raInd] - wi[1]*src[down+step][rbInd] - wi[2]*src[down+step][rcInd] + wi[3]*src[down+step][rdInd];
						a00a = a00a + tmp[0];
						a10a = a10a - tmp[0];
						a01a = a01a + tmp[0];
						a11a = a11a - tmp[0];
						
						tmp[1] = wi[0]*src[down+step][rbInd] + wi[1]*src[down+step][raInd] - wi[2]*src[down+step][rdInd] + wi[3]*src[down+step][rcInd];
						a00b = a00b + tmp[1];
						a10b = a10b - tmp[1];
						a01b = a01b + tmp[1];
						a11b = a11b - tmp[1];
						
						tmp[2] = wi[0]*src[down+step][rcInd] - wi[1]*src[down+step][rdInd] + wi[2]*src[down+step][raInd] - wi[3]*src[down+step][rbInd];
						a00c = a00c + tmp[2];
						a10c = a10c - tmp[2];
						a01c = a01c + tmp[2];
						a11c = a11c - tmp[2];
						
						tmp[3] = wi[0]*src[down+step][rdInd] + wi[1]*src[down+step][rcInd] + wi[2]*src[down+step][rbInd] + wi[3]*src[down+step][raInd];							
						a00d = a00d + tmp[3];
						a10d = a10d - tmp[3];
						a01d = a01d + tmp[3];
						a11d = a11d - tmp[3];
						
						//XXX = XXX.add(t2.multy(a01));
						tmp[0] = wj[0]*src[down][rsaInd] - wj[1]*src[down][rsbInd] - wj[2]*src[down][rscInd] + wj[3]*src[down][rsdInd];
						a00a = a00a + tmp[0];
						a10a = a10a + tmp[0];
						a01a = a01a - tmp[0];
						a11a = a11a - tmp[0];
						
						tmp[1] = wj[0]*src[down][rsbInd] + wj[1]*src[down][rsaInd] - wj[2]*src[down][rsdInd] + wj[3]*src[down][rscInd];
						a00b = a00b + tmp[1];
						a10b = a10b + tmp[1];
						a01b = a01b - tmp[1];
						a11b = a11b - tmp[1];
						
						tmp[2] = wj[0]*src[down][rscInd] - wj[1]*src[down][rsdInd] + wj[2]*src[down][rsaInd] - wj[3]*src[down][rsbInd];
						a00c = a00c + tmp[2];
						a10c = a10c + tmp[2];
						a01c = a01c - tmp[2];
						a11c = a11c - tmp[2];
						
						tmp[3] = wj[0]*src[down][rsdInd] + wj[1]*src[down][rscInd] + wj[2]*src[down][rsbInd] + wj[3]*src[down][rsaInd];	
						a00d = a00d + tmp[3];
						a10d = a10d + tmp[3];
						a01d = a01d - tmp[3];
						a11d = a11d - tmp[3];
						
						//XXX = XXX.add(t1.multy(t2.multy(a11)));
						tmp[0] = wj[0]*src[down+step][rsaInd] - wj[1]*src[down+step][rsbInd] - wj[2]*src[down+step][rscInd] + wj[3]*src[down+step][rsdInd];
						tmp[1] = wj[0]*src[down+step][rsbInd] + wj[1]*src[down+step][rsaInd] - wj[2]*src[down+step][rsdInd] + wj[3]*src[down+step][rscInd];
						tmp[2] = wj[0]*src[down+step][rscInd] - wj[1]*src[down+step][rsdInd] + wj[2]*src[down+step][rsaInd] - wj[3]*src[down+step][rsbInd];
						tmp[3] = wj[0]*src[down+step][rsdInd] + wj[1]*src[down+step][rscInd] + wj[2]*src[down+step][rsbInd] + wj[3]*src[down+step][rsaInd];	
						tmp2[0] = wi[0]*tmp[0] - wi[1]*tmp[1] - wi[2]*tmp[2] + wi[3]*tmp[3];
						tmp2[1] = wi[0]*tmp[1] + wi[1]*tmp[0] - wi[2]*tmp[3] + wi[3]*tmp[2];
						tmp2[2] = wi[0]*tmp[2] - wi[1]*tmp[3] + wi[2]*tmp[0] - wi[3]*tmp[1];
						tmp2[3] = wi[0]*tmp[3] + wi[1]*tmp[2] + wi[2]*tmp[1] + wi[3]*tmp[0];	
						
						a00a = a00a + tmp2[0];
						a00b = a00b + tmp2[1];
						a00c = a00c + tmp2[2];
						a00d = a00d + tmp2[3];
						
						a10a = a10a - tmp2[0];
						a10b = a10b - tmp2[1];
						a10c = a10c - tmp2[2];
						a10d = a10d - tmp2[3];
						
						a01a = a01a - tmp2[0];
						a01b = a01b - tmp2[1];
						a01c = a01c - tmp2[2];
						a01d = a01d - tmp2[3];
						
						a11a = a11a + tmp2[0];
						a11b = a11b + tmp2[1];
						a11c = a11c + tmp2[2];
						a11d = a11d + tmp2[3];
						
						src[down][raInd] = a00a;
						src[down][rsaInd] = a01a;
						src[down+step][raInd] = a10a;
						src[down+step][rsaInd] = a11a;
						
						src[down][rbInd] = a00b;
						src[down][rsbInd] = a01b;
						src[down+step][rbInd] = a10b;
						src[down+step][rsbInd] = a11b;
						
						src[down][rcInd] = a00c;
						src[down][rscInd] = a01c;
						src[down+step][rcInd] = a10c;
						src[down+step][rscInd] = a11c;
						
						src[down][rdInd] = a00d;
						src[down][rsdInd] = a01d;
						src[down+step][rdInd] = a10d;
						src[down+step][rsdInd] = a11d;
					
						//MIRROR
						if (q >= 2){
							int rmirr = mirrorIndex(right, matrixWidth);
							int dmirr = mirrorIndex(down, matrixWidth);
							
							mirror(src, down, right, dmirr, rmirr);
							
							rmirr = mirrorIndex(right + step, matrixWidth);
							mirror(src, down, right + step, dmirr, rmirr);
							
							dmirr = mirrorIndex(down + step, matrixWidth);
							rmirr = mirrorIndex(right, matrixWidth);
							mirror(src, down + step, right, dmirr, rmirr);
							
							rmirr = mirrorIndex(right + step, matrixWidth);
							mirror(src, down + step, right + step, dmirr, rmirr);
						}
						//MIRROR
						
					}
				}
			}
		}
		return src;
	}
	
	protected static void mirror(float[][] src, int down, int right, int dmirr, int rmirr){
		int rmirr4 = 4*rmirr;
		int right4 = 4*right;
		if (rmirr != right){
			src[down][rmirr4] = src[down][right4];
			src[down][rmirr4+1] = src[down][right4+1];
			src[down][rmirr4+2] = -src[down][right4+2];
			src[down][rmirr4+3] = -src[down][right4+3];
		}
		if (dmirr != down){
			src[dmirr][right4] = src[down][right4];
			src[dmirr][right4+1] = -src[down][right4+1];
			src[dmirr][right4+2] = src[down][right4+2];
			src[dmirr][right4+3] = -src[down][right4+3];
		}
		if (rmirr != right && dmirr != down){
			src[dmirr][rmirr4] = src[down][right4];
			src[dmirr][rmirr4+1] = -src[down][right4+1];
			src[dmirr][rmirr4+2] = -src[down][right4+2];
			src[dmirr][rmirr4+3] = src[down][right4+3];
		}
	}
	
	protected static void inverse(float[] arr){
		int dl = arr.length/2;
		int st = arr.length - 1;
		int j = 0;
		for (int i = 0; i < st; i++){
			if (i < j){
				float s0 = arr[i];
				arr[i] = arr[j];
				arr[j] = s0;
			}
			int k = dl;
			while (k <= j){
				j = j - k;
				k = k >> 1;
			}
			j = j + k;
		}
	}
	
	protected static void inverse(float[][] src){
		int N = src.length;
		for (int i = 0; i < N; i++){
			inverse(src[i]);
		}
		float[] arr = new float[N];
		for (int i = 0; i < N; i++){
			for (int j = 0; j < N; j++){
				arr[j] = src[j][i];
			}
			inverse(arr);
			for (int j = 0; j < N; j++){
				src[j][i] = arr[j];
			}
		}
	}
	
	protected static int mirrorIndex(int srcIndex, int matrixWidth){
		int matrixIndex = srcIndex/matrixWidth * matrixWidth;
		return (matrixIndex == srcIndex) ? srcIndex : matrixWidth - srcIndex + 2*matrixIndex;
	}
	
	protected static float[][] transformToGyperComplexMatrix(float[][] src){
		float[][] res = new float[src.length][src[0].length*4];
		for (int i = 0; i < src.length; i++){
			for (int j = 0; j < src.length; j++){
				res[i][j*4] = src[i][j];
			}
		}
		return res;
	}
	
	public static void main(String[] args){
		float[][] arr = Generator.generateSimpleMatrix(2048);
		//Printer.printVector(arr);
		long time = System.nanoTime();
		//Printer.printVector(arr);
		
		
		_fft(arr);
		time = System.nanoTime() - time;
		System.out.println("RESULT: " + time);
		
		
	}
	
}
